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Preface 


The  objective  of  this  new  study  is  an  examination  of  the  influence  of  fracturing  and 
anisotropy  in  the  source  medium  or  near  the  source  on  radiation  patterns  from  explosions, 
including  the  generation  of  S  V  and  SH.  In  previous  work  under  other  contracts  we  have 
studied  radiation  patterns  of  P  and  S  waves  for  explosions  set  off  in  anisotropic  media, 
especially  for  cases  that  mimic  earthquakes  and/or  generate  SH  waves.  This  report  consists 
of  a  manuscript  submitted  to  Geophysical  Journal  International  presenting  a  formalism  for 
calculating  the  effects  of  scattering  from  fracture  zones  of  material,  such  as  might  be  found 
near  the  source  in  a  test  site. 
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ABSTRACT 


We  present  a  ray-Born  method  for  the  computation  of  scattered  wavefields  in  a  general  3- 
D,  ajiisotropic  medium.  This  approach  applies  ray  methods  to  the  computation  of  Green's 
tensors  for  a  background  earth  model  and  uses  the  Born  approximation  to  determine  the 
scattered  wavefield  from  each  volume  element  within  a  discretized  model  of  heterogeneity. 
The  application  of  these  two  approximate  methods  requires  that  the  background  model 
be  relatively  smooth  compared  to  a  wavelength  for  the  validity  of  ray  theory,  but  that 
the  scattering  heterogeneity  have  short  characteristic  length  compared  to  the  propagating 
wavelength  for  accurate  use  of  the  Born  approximation.  Comparisons  of  ray-Born  results  to 
the  complete  solution  for  scattering  from  an  elastic  sphere  show  that  this  method  works  fairly 
well  for  wavelengths  on  the  order  of  five  times  Iwger  than  the  length  scales  typical  of  the 
heterogeneity,  but  then  breaks  down  due  to  the  failure  of  the  Born  approximation.  With  this 
restriction  in  mind,  the  method  is  applied  to  a  hypothetical  layered  earth  model  containing 
a  thin,  laterally  extensive  fracture  zone.  The  results  show  that  scattering  from  shear  waves 
gives  unique  information  on  fracture  orientation,  since  the  properties  of  the  reflected  events 
depend  significantly  on  the  orientation  of  the  incident  S-wave  polarization  with  respect  to  the 
fractures.  For  example,  as  the  polarization  of  an  incident  SH-wave  varies  from  perpendicular 
to  parallel  to  aligned  vertical  fractures,  the  scattered  wave  amplitude  decreases  to  zero.  In 
intermediate  directions,  the  polarization  of  the  scattered  wave  includes  some  SV  energy.  On 
the  other  hand,  compressional  wave  reflections  show  essentially  no  variation  with  orientation 
of  the  fractures.  Modeling  of  waves  scattered  from  fracture  zones  in  V"SP  data  from  the 
Larderello  geothermal  field  in  Italy  demonstrates  the  applicability  of  the  method  to  modeling 
of  field  data  and  suggests  that  at  least  in  this  locality,  anisotropic  fracturing  is  not  completely 
responsible  for  the  observations.  Analysis  of  the  Fresnel  zones  affecting  reflections  from  the 
thin  fracture  zones  responsible  for  the  scattering  allows  a  delineation  of  regions  of  more 
intense  fracturing,  which  is  important  for  the  development  of  geothermal  resources. 
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INTRODUCTION 


The  modeling  of  seismic  wave  propagation  to  understand  the  effects  of  various  earth  struc¬ 
tures  on  observations  and,  conversely,  to  infer  rock  properties  from  data  is  hindered  by  the 
complexity  of  geological  materials.  In  many  cases,  the  geologic  structures  of  interest  are 
small  compared  to  the  three-dimensional  volume  through  which  the  waves  are  transmitted, 
leading  to  practical  difficulties  in  the  implementation  of  numerical  schemes,  such  as  finite 
differences.  Wave  propagation  is  also  complicated  in  many  cases  by  anisotropy,  which  may 
be  an  inherent  anisotropy  of  the  minerals  or  an  effective  anisotropy  due  to  the  pre.scnce  of 
a  stack  of  thin  isotropic  layers  (Crampin,  1981).  Either  of  these  effects  usually  leads  to  a 
transversely  isotropic  medium  with  a  vertical  axis  of  symmetry.  Another  important  source  of 
effective  anisotropy  in  the  crust  is  the  alignment  of  fractures  in  an  otherwise  isotropic  and  ho¬ 
mogeneous  layer,  also  leading  to  a  transversely  isotropic  medium  (Bamford  and  Nunn,  1979; 
Crampin,  1981;  Crampin  et  al.,  1986;  Leary  et  al.,  1987).  Frequently  the  least  principal 
stress  is  known  through  geological,  seismological  or  borehole  observations  to  be  horizontal 
at  depth  (Zoback  and  Zoback,  1980;  Jamison  and  Cook,  1980;  Hickman  et  al.,  1988;  Evans 
et  al.,  1989),  and  this  configuration  results  in  a  vertical  alignment  of  cracks.  The  medium 
then  has  a  horizontal  axis  of  symmetry. 

There  have  been  many  analyses  of  the  effects  of  fractured  layers  of  the  earth  on  synthetic 
seismograms  or  field  data  (Crampin,  1978;  Crampin,  1981;  Crampin  and  McGonigle,  1985; 
Leary  and  Henyey,  1985;  Martin  and  Davis,  1987;  Ben-Menahem  and  Sena,  1990;  Mandal 
and  Toksoz,  1990;  Spencer  and  Chi,  1991;  Mueller,  1991).  Most  of  these  studies  focus  on  the 
effects  of  shear  wave  splitting  and  the  consequent  anomalies  in  the  polarization  of  observed 
shear  wave  data.  These  approaches  are  very  successful  in  locations  where  the  fractured, 
anisotropic  region  is  relatively  thick.  In  these  areas,  the  two  quasi-shear  waves  are  observed 
at  the  receiver  location  with  enough  separation  in  time  that  they  can  be  uniquely  identified. 

However,  it  is  not  realistic  to  suggest  that  the  earth  is  everywhere  uniformly  fractured. 
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Many  practical  problems  require  the  location  of  relatively  small  fracture  zones  within  the 
surrounding  bedrock.  Some  important  examples  of  such  problems  are  the  location  of  frac¬ 
tures  for  the  development  of  geothermal  fields  (Leary  and  Henyey,  1985;  Batini  et  ai.  1985a; 
Batini  et  al.,  1991),  for  nuclear  waste  disposal  (Green  and  Mair,  1983;  Carswell  and  Moon. 
1985),  or  for  other  purposes  (Juhlin  et  al.,  1991).  These  are  examples  of  situations  where 
the  dimensions  of  the  fracture  zones  are  often  not  sufficiently  thick  that  shear  wave  splitting 
will  be  detectable,  even  if  the  zone  is  anisotropic.  Therefore,  a  different  form  of  evidence  for 
the  fracturing  must  be  sought  in  seismic  data. 

To  achieve  this  objective,  we  apply  a  ray-Born  computation  of  synthetic  seismograms 
to  model  the  scattering  of  elastic  waves  by  relatively  small  features  in  three-dimensional 
structures.  This  approach  is  similar  to  that  of  Beydoun  and  Mendes  (1989).  who  outline  a 
ray-Born  algorithm  for  migration  or  inversion  problems  in  three-dimensional,  isotropic  me¬ 
dia  with  examples  of  application  to  two-dimensional  problems.  Two  basic  approximations 
are  employed  to  compute  the  wavefields  expected  for  relatively  complicated  media.  First, 
ray  theory  is  used  to  compute  the  Green’s  tensors  for  a  background  earth  model.  Because 
asymptotic  ray  theory  is  employed,  this  background  reference  model  must  be  sufficiently 
smooth  that  the  high  frequency  assumption  is  valid  (Ben-Menahem  and  Beydoun,  1985). 
Ray  tracing  has  the  advantage  of  allowing  the  computation  of  Green’s  tensors  for  large, 
three-dimensional  models  with  minimal  computation  time.  The  second  step  is  to  use  the 
Born  approximation  to  compute  the  elastic  waves  generated  by  perturbations  to  the  refcr- 
(mce  model,  a  method  which  expresses  the  effects  of  the  perturbations  as  secondary  sources 
radiating  energy  as  they  are  encountered  by  the  wave  propagating  in  the  background  medium 
(Gubernatis  et  ai,  1977;  Wu  and  Aki,  1985;  Gibson  and  Ben-Menahem,  1991).  .Application 
of  this  approximation  is  valid  only  for  small,  short  wavelength  features  of  an  earth  model. 
We  show  how  the  ray-Born  method  can  be  extended  to  fully  general  anisotropic  and  inho¬ 
mogeneous  earth  models  using  the  high  frequency  Green’s  tensor  for  anisotropic  background 
media  (Ben-Menahem  et  al.,  1990)  and  the  expressions  for  the  secondary  source  radiation 


by  anisotropic  perturbations  developed  in  Gibson  and  Ben-Menahem  (1991).  We  apply  the 
method  to  three-dimensional  isotropic  earth  models  with  anisotropic  perturbations. 

.\fter  describing  the  ray-Born  algorithm,  including  a  brief  review  of  the  Born  approxima¬ 
tion,  we  explore  the  question  of  the  accuracy  of  the  Born  approximation  in  computing  the 
scattered  waves  caused  by  varying  degrees  of  heterogeneity.  This  is  accomplished  by  com¬ 
paring  the  ray-Born  results  with  the  complete  solution  for  waves  scattered  from  an  elastic 
sphere  by  an  incident  plane  P-wave.  Using  the  guidelines  for  application  of  the  ray-Born 
algorithm  developed  by  this  comparison,  we  next  examine  the  wavefields  generated  by  a 
thin  but  laterally  extensive  fracture  zone  in  a  hypothetical  layered  earth  model  to  gain  an 
understanding  of  the  effects  of  fracture  orientation  on  seismic  waves  in  an  ideal  case.  Lastly, 
we  examine  a  set  of  VSP  data  from  the  Larderello  geothermal  field  in  Italy,  where  an  impor¬ 
tant  objective  in  the  course  of  exploitation  of  geothermal  resources  is  the  delineation  of  thin 
fracture  zones  several  kilometers  in  depth.  We  present  a  fracture  zone  model  that  explains 
the  observed  data  and  discuss  the  implications  of  the  results  to  determine  fracture  alignment 
in  the  Larderello  field  bcised  on  the  synthetic  results  from  our  model. 

COMPUTATION  OF  THE  BORN  SCATTERED  FIELD 
The  Born  Integral  Equation 

The  vector  form  of  the  source  free  equation  of  motion  for  an  inhomogeneous  medium  is 

^u-V-T  =  0,  (1) 

where  p  is  density,  u  is  the  displacement  vector,  and  T  is  the  stress  tensor  (Ben-Menahem 
and  Singh,  1981).  Derivatives  with  respect  to  time  are  indicated  by  the  dot  symbols  over 
the  vector  u.  Hooke’s  law  relates  the  stress  tensor  and  the  strain  tensor  E  through  the 
fourth-order  elastic  tensor  c: 

T  =  c  :  E.  (2) 
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The  elastic  tensor  obeys  the  usual  symmetry  relationships 


^ijkl  —  ^kttj  —  ^jikl  —  ^ijlki  (3) 

resulting  in  a  total  of  21  independent  elastic  constants. 

The  Born  approximation  to  the  scattered  field  generated  by  heterogeneous  materials  is 
obtained  by  considering  perturbations  of  the  elastic  properties  of  a  prescribed  background 
reference  model: 

Mx)  =  /)°(x)  +  <5p(x) 

c.jfc/(x)  =  c°jfc,(x)  +  ^Cijfc,(x).  (4) 

Here  the  superscript  ®  indicates  a  value  of  the  known  background  model  and  x  is  the  position 
vector.  The  values  6p(x)  and  8c,jki(x)  are  the  perturbations  to  the  reference  values  /)°(x) 
and  c°jj.^(x).  Both  the  background  and  perturbation  models  are  taken  to  be  functions  of 
location  x.  Derivation  of  expressions  for  the  scattered  field  proceeds  by  assuming  that  the 
displacement  field  in  the  total  medium  is  also  given  by  the  sum  of  the  field  u®(x)  which 
would  propagate  in  the  background  material  and  a  scattered  field  <5u(x): 

u(x,t)  =  u°(x,t)  +  6u(x,t).  (-5) 

Gibson  and  Ben-Menahem  (1991)  show  that  substitution  of  equations  (4)  and  (5)  into  the 
equation  of  motion  (1),  neglecting  second-order  terms,  yields  a  solution  for  the  scattered 
field 

6u{x,t)  =  j  dV{x')  J  dt'  G^{x,t:x\t')  ■  {dpu°(x',t')) 

+  J  dVix')  J  dt’  [6c{x',t’)  :  E°(x',t')]  :  VGO{x,  f;  x',  t').  (6) 

The  ~  symbol  indicates  transpose  of  the  Green’s  tensor  G®.  This  Born  approximation  pro¬ 
vides  an  approximation  to  the  scattered  field  in  terms  of  the  incident  wavefield  u®(x.  t)  prop¬ 
agating  in  the  background  medium,  as  E°(x',  t')  is  the  strain  associated  with  this  displace¬ 
ment.  G”(x,  <;  x',  t')  is  the  Green’s  tensor  for  the  background  earth  model,  with  component 


G’°  (x.  representing  the  i  component  disturbance  at  location  x  from  a  source  applied 

at  x'  in  the  xj  direction.  The  perturbations  to  the  material  quantities  act  as  secondary 
sources  for  displacement  fields  propagating  in  the  reference  model,  where  the  product  of 
density  perturbation  and  acceleration  vector  yields  a  single-force  source  vector  at  each  point 
x'  and  time  t'.  Likewise,  perturbations  to  the  elastic  tensor  c  result  in  double-force  source 
terms  which  are  doubly  contracted  against  G”(x,Z;x',  <').  This  moment  tensor  type  source 
=  ScijkiEici  is  more  easily  interpreted  by  writing  it  in  the  following  form,  employing  the 
standard  Voigt  notation: 


'  «5A/„  ■ 

SCii  SC12  SCi,\  SO\5 

6M22 

6C22  ^(^23  ^(^24  <5G25  <5(726 
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In  this  expression,  e°  are  the  components  of  the  strain  tensor  E°(x',t'). 
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Evaluation  of  the  Integral 

Using  equation  (6),  the  Born  scattered  field  for  prescribed  perturbations  6p  and  Scijkt  can 
be  computed  once  the  incident  field  u”(x,<)  and  its  associated  strain  E°(x',t')  are  known. 
The  Green’s  tensor  G°(x,t;x',  t')  is  the  fundamental  quantity  which  must  be  obtained,  as  it 
can  be  used  to  compute  these  background  fields  as  well  as  the  displacements  generated  by 
the  secondary  sources  on  the  right-hand  side  of  equation  (6).  Therefore,  application  of  the 
integral  solution  first  requires  the  knowledge  of  the  Green’s  tensors  for  propagation  in  the 
background  model  and  then  an  algorithm  for  the  actual  evaluation  of  the  integral. 

Following  Beydoun  and  Mendes  (1989),  we  proceed  by  discretizing  the  perturbed  volume 
of  the  earth  model  (Figure  1).  Each  unit  cell  within  the  discretization  is  a  rectangular 
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prism  with  dimensions  dxj,  dx^,  and  dx^.  As  long  as  the  wavelength  of  the  incident  wave  is 
much  greater  than  the  largest  dimension  di,-,  the  scattering  due  to  each  individual  elemental 
volume  is  equivalent  to  Rayleigh  scattering,  which  reduces  the  integration  over  volume  in 
equation  (6)  to  a  simple  multiplication  by  the  elemental  volume  dV  =  dxxdx^dx^  for  a  given 
unit  cell: 

(5u(x,  i)  =  <51/  J  dt'  G°(x,/;x', /')  •  ((5/Ju°(x',  <')) 

+  6V  J  dt'  (6c(0,  t')  :  E°(0,  t'))  :  VG®(x,  <;  0,  t'),  (8) 

or,  in  component  notation, 

Sui(x,t)  =  6V  j  dt'  G,j{x,t;x',t'){6puj^{x',t')) 

+  6VJ  dt'  (,5c,fc,™(x^Oc?„(x^O)G®  .*(x,<;x^<')•  (9) 

Under  the  condition  of  Rayleigh  scattering,  each  cell  in  the  perturbed  volume  therefore  acts 
as  a  point  source  located  at  the  center  of  the  cell,  indicated  by  the  dots  in  Figure  1.  The  point 
source  at  each  lattice  point  is  equivalent  to  some  combination  of  single  or  double  forces  as 
presented  in  equation  (6).  Evaluation  of  the  integral  is  most  easily  accomplished  by  simply 
summing  the  contributions  of  each  of  the  point  sources  within  the  scattering  volume,  which 
corresponds  to  the  most  elementary  implementation  of  the  definition  of  an  integral  in  terms 
of  the  limit  of  a  sum  of  increasingly  finely  discretized  representations  of  the  range  of  the 
integral  (Budak  and  Fomin,  1983).  As  long  as  the  discretization  of  the  volume  is  sufficiently 
small,  the  evaluation  of  the  integral  by  this  means  should  be  sufficiently  accurate. 

A  constraint  on  the  sizes  of  the  discretization  intervals  dxi,  dx2  and  dx^  is  made  by 
considering  the  dominant  period  of  the  incident  wave  and  therefore  the  incident  wavelength. 
Adequate  construction  of  the  scattered  wave  requires  at  least  four  samples  per  period  of  the 
wave.  For  a  reflected  (back  scattered)  wave,  the  interval  between  the  signals  by  adjacent 
pixels  will  be  2dx,/u  for  a  wave  propagating  in  the  x,  direction  at  velocity  r.  To  obtain  an 
adequate  sampling  then  requires  that  the  spa.cing  dx,  be  less  than  or  equal  to  1/8  of  the 
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wavelength.  For  nearly  vertically  travelling  waves,  this  constraint  will  apply  to  dxs.  The 
intervals  dxi  and  dx2  can  be  somewhat  larger,  as  the  wavefront  will  be  essentially  tangent 
to  the  X1-X2  plane,  and  the  integrand  in  equation  (6)  will  vary  much  more  slowly  in  these 
directions. 

The  application  of  this  integration  scheme  then  requires  a  knowledge  of  the  Green’s 
tensors  corresponding  to  the  wave  propagation  from  the  primary  source  to  each  elemental 
scattering  volume,  yielding  the  properties  of  each  secondary  point  source,  and  also  from 
the  scattering  volume  to  the  receivers,  which  determines  the  scattered  field  from  each  point 
source.  Ray  methods  provide  a  fast  and  flexible  means  for  performing  these  calculations  in 
general  three-dimensional  layered  models.  The  principal  requirement  for  ray  solutions  to  be 
applicable  is  that  the  wavelength  must  be  much  less  than  the  characteristic  length  scale  of 
the  background  earth  model  (Ben-Menahem  and  Beydoun,  1985).  For  inhomogeneous  and 
anisotropic  media  meeting  this  length  scale  requirement,  the  ray  theoretical  Green’s  tensor 
is  given  by  (Ben-Menahem  et  ai,  1991) 


[g(x)g(xo)]6(T(x|xo)),  (10) 


G°(x,t;x',F)  =  r  w..  [g(x)g(xo)]6(T(x|xo)),  (10) 

4ff[i;(xo)]'  /j(x)i;(x)J(x) 

where  Xo  is  the  source  position,  x  is  the  observation  point,  r;(x)  is  phase  velocity,  and 
u(x)J(x)  is  the  Jacobian  of  the  transformation  from  Cartesian  coordinates  i,  to  ray  coordi¬ 
nates  7j,  where  J(x)  is  given  by 

J(x)  =  —  X  —  .  (11 

dii  572 


The  ray  coordinates  are  specified  to  be  the  two  take-off  angles  of  the  ray  at  the  source,  the 
declination  angle  and  the  azimuthal  angle  (j)  (Figure  2).  It  is  important  to  realize  that 
this  Green’s  tensor  contains  only  quantities  obtained  in  the  course  of  normal  ray  tracing 
algorithms.  It  contains  the  effects  of  geometrical  spreading  on  amplitudes  through  the  ratio  of 
Jacobians  r(xo)o/(xo)/t’(x)7(x),  cis  well  cis  the  material  properties  at  the  source  and  receiver 
points.  The  scalar  amplitude  multiplies  a  dyad  given  by  the  outer  product  of  the  polarization 
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vectors  at  the  source  and  at  the  observation  position,  reproducing  the  usual  reciprocity  of 
elastic  wave  propagation  in  inhomogeneous  media,  whereby  exchange  of  source  and  recei\'er 
positions  result  in  the  equivalence  (Ben-Menahem  and  Singh,  1981) 


(12) 


This  ray  theoretical  Green’s  tensor  will  not  be  able  to  model  aspects  of  wave  propagation  such 
as  caustics  or  shadow  zones,  a  fundamental  limitation  of  the  high  frequency  approximation 
(Cerveny,  1985). 

Given  the  general  ray  theoretical  Green’s  tensor  in  equation  (10),  the  ray- Born  method 
could  be  applied  to  general  anisotropic  media  with  anisotropic  inhomogeneity.  However,  as  a 
first  step,  we  consider  only  an  isotropic,  inhomogeneous  background  model  with  anisotropic 
inclusions  in  order  to  develop  an  understanding  of  the  effects  of  localized  anisotropic  regions 
on  elastic  wave  propagation.  This  also  allows  the  utilization  of  the  elegant  dynamic  ray 
tracing  (DRT)  techniques  in  the  ray-ccntered  coordinates  qi  with  basis  vectors  e/  (Figure  2) 
(Cerveny,  1985).  The  DRT  involves,  in  addition  to  the  standard  determination  of  ray  path 
and  travel  time  (kinematic  ray  tracing,  KRT),  the  integration  of  eight  additional  ordinary 
differential  equations  to  obtain  the  2  by  2  matrices  Q  and  P.  These  components  of  these 
matrices  can  be  expressed  as 


QlJ 

Pu 


dqi 
d-yj 
1  dQjj 
V*  dr 


(13) 

(14) 


Here  r  is  the  travel  time  along  the  ray,  and  l.J  =  1,2.  The  Pjj  are  needed  only  to  obtain 


QiJ- 

The  matrix  Q  is  related  to  the  curvature  of  the  wavefront  (Cerveny,  1985).  Knowledge  of 
the  components  of  this  matrix  yields  geometric  spreading,  which  is  proportional  to  det  Q.  Q 
is  also  used  for  the  application  of  the  paraxial  method,  which  allows  extrapolation  of  travel 
time  and  polarization  vectors  from  a  central  ray  which  has  been  obtained  by  integration 
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to  nearby  observation  points.  This  results  in  significant  savings  in  computation  time,  as 
the  two-point  ray  tracing  is  avoided.  In  an  application  such  as  the  ray-Born  method  where 
the  incident  wavefield  must  be  known  at  a  large  number  of  points,  this  is  an  especially 
valuable  feature.  These  paraxial  ray  tracing  procedures  are  discussed  in  many  references 
(e.g.,  Cerveny,  1985). 

In  addition  to  these  more  typical  paraxial  corrections,  we  also  include  a  correction  to  the 
geometrical  spreading  amplitude  factor  which  can  be  derived  using  geometrical  arguments. 
The  integration  of  both  the  KRT  and  DRT  equations  is  dependent  on  the  choice  of  initial 
values,  but  once  these  values  are  selected  the  integration  of  the  ray  equations  governs  all 
propagation  effects  along  the  ray  path  from  the  initial  point,  within  the  limitations  of  the  va¬ 
lidity  of  ray  methods.  Because  the  initial  conditions  can  be  chosen  along  an  initial  wavefront 
zis  easily  as  from  a  source  point,  it  follows  that  the  additional  geometrical  spreading  from  the 
wavefront  at  central  ray  point  x'  and  time  t(x')  to  the  wavefront  at  observation  point  x  and 
time  t(x)  -i-  Ar  is  in  a  homogeneous,  isotropic  material  equivalent  to  the  distance  between 
wavefronts  uAr.  The  total  geometrical  spreading  is  then 

detQ-fwAT.  (15) 

In  an  inhomogeneous  material,  this  correction  will  be  less  accurate,  but  since  the  paraxial 
corrections  can  only  be  applied  a  relatively  short  distance  from  the  central  ray,  it  should 
not  be  a  severe  limitation;  application  of  this  correction  can  be  restricted  to  homogeneous 
regions  of  a  model.  It  should  also  be  noted  that  in  this  formulation,  the  changes  in  reflection 
and  transmission  coefficients  with  the  variation  in  incidence  angles  on  interfaces  from  the 
central  ray  to  the  ray  which  would  actually  intersect  the  observation  point  are  not  included. 
This  also  should  not  be  a  severe  limitation  except  possibly  near  critical  angles. 

.'\nother  useful  aspect  of  the  DRT  is  that  the  partial  derivatives  may  be  used  in  an 
iterative  two-point  ray  tracing  scheme  for  rays  propagating  from  the  source  to  fixed  receiver 
locations.  While  this  is  not  necessary  for  the  calculation  of  synthetic  seismograms,  it  is  often 
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desirable  for  the  determination  of  accurate  ray  paths  and  will  yield  more  accurate  results 
than  the  straight  application  of  the  paraxial  method.  Once  a  ray  is  traced  with  an  initial 
source  point  x  and  take-off  angles  7°  and  7°,  the  71  and  72  corresponding  to  the  ray  arriving 
at  the  desired  observation  point  can  be  estimated  through 

11^1°  +  (Q~^)ijqj-  (16) 

We  apply  this  result  to  two  steps  in  the  ray-Born  procedure.  First,  for  the  computation  of  the 
background  synthetic  seismograms,  we  perform  an  iterative,  two- point  ray  tracing  procedure 
where  we  shoot  a  fan  of  rays  over  some  range  of  prescribed  take-off  angles.  Beginning  with 
the  closest  ray  to  each  receiver,  we  repeatedly  apply  equation  ( 16)  until  the  ray  arrives  with 
some  distance  of  the  desired  receiver  point.  For  layered  models,  the  ray  arriving  essentially 
at  the  receiver  point  can  usually  be  determined  in  three  or  fewer  iterations,  providing  a 
very  rapid  determination  of  ray  paths  for  multiply  reflected  and  transmitted  rays.  For  more 
complicated  models,  the  procedure  may  not  converge  well,  in  which  case  the  straight  paraxial 
method  can  be  applied.  The  second  application  of  equation  (16)  is  in  the  computation  of 
the  Green’s  tensor.  Since  the  dyad  in  equation  (10)  contains  the  polarization  vector  at  boih 
the  source  and  receiver  points,  to  compute  the  tensor  using  the  paraxial  method  at  the 
scattering  lattice  we  employ  (16)  to  correct  the  polarization  at  the  start  of  the  central  ray 
to  be  approximately  that  of  the  ray  joining  the  source  point  and  scattering  point.  Cerveny 
et  al.  (1987)  outline  the  conversion  of  the  results  for  the  polarization  dyad  in  ray-centered 
coordinates  to  Cartesian  coordinates.  Both  SH  and  SV  waves  are  automatically  included  in 
this  procedure  applied  in  ray-centered  coordinates. 

The  algorithm  can  be  summarized  as  follows: 

I 

I 

1.  Synthetic  seismograms  are  computed  for  the  background  model  using  the  iterative 
two-point  ray  tracing  technique  or  the  paraxial  method.  Green's  tensors  are  com¬ 
puted  and  then  the  desired  primary  source  type  is  applied  to  compute  the  background 
displacement  field. 
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2.  Green’s  tensors  are  conij)uted  for  pro[)agation  from  the  primary  source  to  eacfi  ele¬ 
mental  volume  in  the  scattering  region.  The  Green’s  tensors  are  also  computed  for 
propagation  from  the  scattering  lattice  to  receivers.  It  is  generally  more  convenient  to 
do  the  ray  tracing  from  receivers  to  the  scattering  volume  and  then  apply  the  princi¬ 
ple  of  reciprocity,  transposing  the  Green’s  tensors  (Ben-Menahem  and  Singh,  1981). 
Paraxial  corrections  are  applied  to  rays  passing  near  scattering  points  to  minimize  the 
amount  of  ray  tracing  necessary. 

3.  Perturbations  to  the  material  properties  are  specified,  and  the  integration  of  equation 
(6)  is  performed  by  summing  the  contribution  of  each  point  source  within  the  scattering 
lattice.  The  scattered  field  resulting  from  this  calculation  is  then  added  to  the  primary 
field  to  produce  the  final  synthetic  seismograms. 

An  advantage  of  this  method  is  that  the  Green’s  tensors  may  be  saved  prior  to  the  application 
of  the  primary  source  of  the  Born  approximation.  This  makes  possible  rapid  comparison  of 
the  effects  of  different  primary  sources  or  different  material  perturbations  on  the  wavefields, 
as  the  ray  tracing,  the  most  time  consuming  part  of  the  algorithm,  need  only  be  done  once. 

COMPARISON  OF  RAY-BORN  RESULTS  TO  A  KNOWN 

ANALYTIC  SOLUTION 

The  preceding  ray- Bom  algorithm  is  in  principle  very  general  and  can  be  applied  to  a  wide 
variety  of  three-dimensional  earth  models,  both  isotropic  and  anisotropic.  As  long  as  the 
ray  tracing  can  be  satisfactorily  accomplished,  the  integral  over  volume  of  heterogeneity  is 
simply  computed  by  summation.  However,  prior  to  application  of  the  method  to  general 
problems,  it  is  desirable  to  have  some  knowledge  ais  to  the  accuracy  and  validity  of  the 
resulting  synthetic  seismograms.  For  this  purpose,  we  compare  the  ray-Born  scattered  field 
results  for  a  spherical  object  to  the  known  solution  in  terms  of  spherical  harmonics.  This 
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comparison  was  chosen  because  it  provides  a  true  3-D  test  of  the  ray-Born  method  for  an 
obstacle  of  finite  extent.  If  the  heterogeneity  were  of  infinite  volume,  such  as  a  thin  layer, 
then  a  synthesis  of  ray-Born  results  in  practical  implementations  would  require  a  truncation 
of  the  range  of  the  integral  in  equation  (6)  to  some  finite  portion  of  the  thin  layer.  In  this 
case,  the  comparisons  would  have  to  simultaneously  examine  both  the  accuracy  of  the  ray- 
Born  method  itself  and  the  validity  of  truncating  the  range  of  integration,  complicating  the 
conclusions  that  could  be  drawn  from  the  results. 

The  scattered  displacement  field  generated  by  a  monochromatic  plane  wave  incident  on  an 
isotropic,  elastic  sphere  in  an  infinite  homogeneous  medium  was  obtained  by  Ying  and  Truell 
(1956).  Derivation  of  this  solution  begins  by  expressing  the  components  of  displacement  for 
a  planar  P-wave  incident  along  the  z  axis, 

V.(i3,<)  = 

V,(X3,0  =  (17) 

in  terms  of  an  expansion  of  a  potential  in  spherical  harmonics: 

Ui  =  -V0, 

1 

0,  =  —  -I- l)j„,(^-,r)F„(cos^).  (18) 

m=0 

In  these  expressions,  ki  =  wq,  where  q  is  the  compressional  wave  velocity  in  the  infinite 
medium,  jm  is  the  spherical  Bessel  function  of  the  first  kind,  and  is  the  Legendre  polyno¬ 
mial.  Spherical  coordinates  r  and  0  are  illustrated  in  Figure  3.  The  scattered  displacements 
are  given  by 

V,  =  v,e'"‘ 

V,  =  -V0, -f  V  X  V  X  (rrll,),  (19) 

where  i/>j  and  II,,  are  the  compressional  and  shear  wave  potentials  respectively,  with  harmonic 
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expansions 


oo 

0,  =  ^  Amhmikir)Pmicos6) 

m=0 

oo 

n.  =  ^  B^/i„(/c.r)P„(cos0).  (20) 

m=0 

Here  hm  is  the  spherical  Bessel  function  of  the  third  kind,  Kj  =  u}0  is  the  shear  wave  number, 
and  Am  and  Bm  are  the  as  yet  unknown  coefficients  of  the  expansion.  The  displacement  field 
inside  the  sphere  is  expressed  in  similar  expansions,  giving  two  more  unknown  coefficients 
Cm  and  Dm-  Boundary  conditions  on  continuity  of  stress  and  displacement  are  applied  at 
the  surface  of  the  spherical  obstacle,  resulting  in  four  simultaneous  equations  that  must  be 
solved  for  the  coefficients  in  each  term  m  in  the  expansions. 

Since  a  general  seismic  signal  is  actually  composed  of  contributions  from  multiple  frequen¬ 
cies,  we  apply  a  discrete  wavenumber  algorithm  to  compute  synthetic  seismograms  as  follows. 
A  range  of  discrete  frequencies  is  specified,  fi  =  A/,2A/,  3A/, ..., /nyq,  where  /„y,  =  1/2A< 
is  the  Nyquist  frequency  corresponding  to  the  time  domain  sample  interval  At.  Each  fre¬ 
quency  fi  corresponds  to  a  wavenumber  The  boundary  value  problem  is  then  solved  at 
frequencies  fi  for  the  coefficients  Am,  Bm,  Cm  and  Dm-  At  each  frequency,  the  harmonic 
expansion  is  carried  out  to  sufficiently  high-order  m„  that  coefficients  Am„,  Bm„,  Cm„  and 
Dm„  are  negligible  compared  to  Aq,  Bq,  Co  and  Dq.  The  Bessel  functions  of  the  first  and  sec¬ 
ond  kinds  of  arbitrary  order  are  computed  using  a  Miller’s  algorithm  appropriately  modified 
for  spherical  Bessel  functions  (Press  et  ai,  1988).  Finally,  after  the  scattered  field  r  and  0 
components  have  been  computed  at  each  frequency,  this  impulse  response  is  convolved  with 
the  spectrum  of  the  incident  plane  wave  signal  through  frequency  domain  multiplication, 
and  an  inverse  FFT  is  applied  to  produce  the  time  domain  response.  In  this  way,  the  full 
waveform  elastic  wave  scattered  field  of  the  spherical  obstacle  can  be  computed  for  a  general 
incident  plane  wave  pulse. 

comparison  of  the  ray-Born  and  discrete  wavenumber  algorithm  results  was  made  using 
a  spherical  object  of  radius  a  =  0.5  km  surrounded  by  a  ring  of  receivers  at  a  distance  of 
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40  km  (Figure  4).  This  configuration  of  receivers  allows  a  comparison  of  both  forward  (90°) 
and  back  scattering  results  (270°).  The  infinite  medium  velocities  were  Qq  =  4.5  km/sec  and 
0Q  =  2.7  km/sec,  with  density  2.67  g/cm^;  velocities  within  the  spherical  region  were  set 
to  a  =  4.545  km/s  and  0  =  2.727  km/s,  perturbations  of  1%.  Density  was  kept  constant, 
so  that  only  the  elastic  Lame  parameters  were  varied.  For  the  ray-Born  calculations,  a 
10  X  10  X  20  lattice  with  spacings  dx\  =  dxi  =  0.10  km  and  di-^  =  0.05  km  was  specified 
and  centered  on  the  origin  (note  that  ij,  and  in  correspond  to  i,  y  and  z  coordinates  in 
Figure  4).  A  finer  spacing  in  the  2  direction  was  used  for  accurate  evaluation  of  the  ray-Born 
integral  since  this  was  the  direction  of  the  incident  wave  propagation  and  thus  the  most  rapid 
variation  of  the  integrand.  Amplitude  and  phase  of  the  incident  plane  wave  are  constant  in 
both  X  and  y.  The  perturbations  at  all  lattice  points  at  a  distance  larger  than  0.5  km  were 
set  to  zero  providing  an  approximation  to  the  spherical  in  homogeneity.  For  both  methods, 
the  source  wavelet  used  was  given  by 

s(L)  =  cos(wo(t  -  to)),  (21) 

where  to  is  the  arrival  time,  a;o  =  27r/o  contains  the  center  frequency  /o,  and  7  is  a  free 
parameter  which  we  set  to  three.  The  resulting  radial  component  P-wave  synthetic  seismo¬ 
grams  from  the  two  methods  for  center  frequencies  of  0.25,  1,  2.5  and  5  Hz  are  compared  in 
Figures  5  and  6.  The  important  parameter  to  determine  the  validity  of  the  Born  approxi¬ 
mation  is  the  ratio  of  wavelength  to  sphere  diameter  77  =  A./d.  These  frequencies  we  have 
examined,  0.25,  1,  2.5  and  5  Hz,  provide  ratios  rjp  =  18.  4.5,  1.8  and  0.9.  respectively,  for 
compressional  waves.  Examination  of  the  results  for  ijp  =  ]S  (Figures  5 A.  6A)  shows  that 
the  two  solutions  are  very  similar,  which  is  not  surprising  as  this  long  wavelength  reproduces 
the  Rayleigh  scattering  result,  where  the  forward  (90°)  and  back  scattering  (270°)  ampli¬ 
tudes  are  the  same.  A  plot  of  the  maximum  amplitudes  as  a  function  of  scattering  direction 
shows  that  there  is  a  systematic  difference  in  the  two  results,  with  the  ray-Born  amplitude 
about  15%  to  20%  less  than  the  di.scretc  wavenumber  result  (Figure  7).  An  error  of  ap- 
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proximately  2%  results  simply  because  the  volume  of  the  discretized  version  of  the  sphere  is 
98%  of  the  true  sphere.  The  rest  of  the  variation  is  explained  by  complicated  interactions 
of  the  plane  wave  with  the  sphere.  Because  the  impulse  response  of  the  sphere  increases 
more  rapidly  with  frequency  than  the  amplitude  of  the  incident  wavelet  decreases  for  the 
low  frequency  source  wavelet,  the  maximum  response  after  convolution  is  actually  shifted  to 
a  slightly  higher  frequency  (Gibson,  1991).  As  the  amplitude  of  Rayleigh  scattering  goes  as 

this  increases  the  amplitude  of  the  discrete  wavenumber  solution.  This  is  not  reflected 
in  the  ray- Born  solution,  since  it  only  includes  the  u,’o  specified  for  the  center  frequency  of 
the  incident  wavelet. 

As  the  incident  center  frequency  increases  to  1  Hz  and  the  ratio  rip  decreases  to  4.-5  (Fig¬ 
ures  5B,  6B).  the  two  methods  compare  about  cis  well  as  for  7/p=18,  noting  also  the  amplitude 
plot  in  Figure  7.  Errors  are  less  than  20%  in  all  directions  for  the  P-wave  amplitudes  and  less 
than  10%  for  all  back  scattered  energy.  Decreasing  the  incident  wavelength  further,  however, 
causes  the  results  for  the  ray- Born  method  to  begin  to  degrade.  We  see  that  the  spherical 
harmonic  solution  predicts  that  the  back  scattering  amplitude  at  270®  will  become  approx¬ 
imately  constant  for  shorter  wavelengths  although  forward  scattering  amplitude  incre^lses 
significantly  with  frequency  (Figures  5C,  6C).  In  one  sense,  the  comparison  of  maximum 
amplitudes  is  not  a  strictly  valid  measure  of  equivalence  for  the  higher  frequency  results 
in  back  scattering,  since  the  discrete  wavenumber  results  show  that  for  rjp  —  0.9  (5  Hz), 
both  shear  and  compressional  wave  back  scattering  have  the  form  of  two  reflections  from  the 
front  and  back  of  the  sphere,  while  the  ray-Born  method  yields  only  a  single  wide  pulse  due 
to  insufficient  cancellation  of  back  scattered  waves.  Therefore,  comparing  the  amplitude  of 
the  single  ray-Born  wavelet  with  the  values  for  the  two  reflected  arrivals  from  the  discrete 
wavenumber  results  does  not  reflect  the  true  mismatch  of  waveforms.  The  scattered  shear 
waves  were  shown  by  Gibson  (1991)  to  lead  to  the  same  conclusions  regarding  the  accuracy 
of  the  ray-Born  method. 

It  is  worth  noting,  however,  that  in  spite  of  the  limitations  of  the  Born  approximation 


the  results  for  forward  scattering  at  90°  from  both  solution  methods  still  tend  to  compare 
relatively  well  at  higher  frequencies.  Also  encouraging  is  that  both  methods  predict  a  much 
shorter  pulse  of  higher  apparent  frequency  for  forward  scattering  than  for  back  scattering  and 
that  the  total  width  in  time  of  energy  at  each  observation  point  is  the  same.  The  ray- Born 
method  succeeds  in  matching  the  gross  features  of  the  wavefields,  including  achieving  at  least 
some  of  the  general  trends  of  amplitude  variation  with  respect  to  both  scattering  direction 
and  increasing  frequency.  Details  of  the  reflections  are  missing  from  the  Born  approximation. 

We  computed  the  scattered  fields  using  both  methods  for  P  and  S-wave  velocity  pertur¬ 
bations  of  10%  and  50%  as  well  as  the  1%  perturbation  results  shown  here.  The  comparisons 
are  essentially  the  same  as  these  results,  although  the  details  of  the  complete  waveform  so¬ 
lution  change  slightly  for  the  50%  velocity  perturbation.  These  are  not  reproduced  by  the 
ray-Born  method,  which  has  strictly  linear  behavior  with  respect  to  variation  in  elastic  con¬ 
stants.  The  amplitude  comparisons  become  only  slightly  worse  as  the  velocity  perturbations 
increase  to  50%.  'I'lic  more  significant  failures  of  the  ray-Born  solution  are  still  caused  by 
increasing  frequency.  Therefore,  we  conclude  that  a  principal  requirement  for  accurate  solu¬ 
tions  using  the  ray-Born  technique  with  the  degrees  of  heterogeneity  we  considered  is  that 
the  incident  and  scattered  wavelengths  must  be  four  to  five  times  larger  than  the  length  scale 
of  the  heterogeneity,  a  restriction  which  for  the  case  above  corresponds  to  frequencies  less 
than  or  equal  to  1  Hz.  As  the  magnitude  of  the  perturbations  decreases,  this  requirement 
can  be  relaxed  as  the  Born  approximation  will  be  more  accurate  for  weaker  heterogeneity. 
Conversely,  for  equal  accuracy  for  larger  perturbations,  a  somewhat  larger  ratio  of  wave¬ 
length  to  heterogeneity  dimension  must  be  imposed  as  the  stronger  velocity  changes  reduce 
the  accuracy  of  the  Born  approximation.  Even  when  the  wavelength  ratio  is  larger  than  5. 
subtle  aspects  of  wave  interaction  with  heterogeneity  can  degrade  results,  as  for  the  0.25  Hz 
scattering  for  our  spherical  model,  where  the  resonances  of  the  sphere  cause  a  shift  of  the 
principal  response  frequency. 

The.se  conclusions  regarding  the  accuracy  and  validity  of  the  ray-Born  method  are  sup- 


ported  by  the  analysis  of  Hudson  and  Heritage  (1981).  By  assuming  that  the  quantities 
|<5A/A“|,  \Spl p°\,  1^x1,  and  are  all  small,  it  was  shown  that  the  series  solution 

returns  the  Born  approximation.  Therefore,  conditions  for  validity  of  the  Born  approxima¬ 
tion  according  to  Hudson  and  Heritage  (1981)  can  be  stated  as: 
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Here  A^  is  shear  wavelength  in  the  background  medium.  These  conditions  are  very  similar  to 
those  discussed  above  regarding  results  from  the  comparisons  of  the  more  complete  numerical 
solutions  with  ray- Born  results. 

Although  density  perturbations  were  not  considered  in  the  above  examples  since  our  ap¬ 
plications  consider  only  elastic  constant  perturbations,  the  restrictions  on  accuracy  of  the 
ray- Born  method  should  be  much  the  same  as  for  elastic  properties.  Both  types  of  perturba¬ 
tions  occur  in  convolutional  integrals  over  volume  and  time  of  the  same  form  in  equation  (6). 
Therefore,  the  main  factor  that  controls  the  accuracy  of  the  Born  approximation  is  still  the 
rapidity  of  variation  of  the  incident  wavefield  with  respect  to  the  volume  of  the  heterogeneity. 
The  only  difference  is  that  for  density,  the  particle  acceleration  multiplies  the  perturbation, 
whereeis  the  strain  multiplies  the  elcistic  constant  variation.  Both  acceleration  and  strain 
vary  with  the  same  incident  wavelength. 


S  i  NTHETIC  SEISMOGRAMS  FROM  A  FRACTURE  ZONE 

MODEL 


Given  the  guidelines  established  in  the  comparisons  with  the  complete  solution  for  scattering 
by  an  elaistic  sphere,  we  apply  the  method  to  a  layered  earth  model  containing  a  relatively 
thin  but  laterally  extended  zone  of  fracturing  in  the  subsurface.  Using  the  ray-Born  method, 
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we  can  easily  compare  the  expected  seismic  response  of  the  model  without  the  fractured 
region  and  with  isotropic  and  anisotropic  fracturing  present.  In  addition,  the  effects  of 
various  crack  filling  materials  on  the  back  scattered  displacement  fields  are  tested. 

The  background  earth  model  is  presented  in  Figure  8.  It  consists  of  three  layers  overlying 
a  half-space.  The  source  is  located  at  a  depth  of  0.25  km  in  the  top  layer,  which  is  0.50  km 
thick,  and  receivers  are  located  on  the  free  surface  at  offsets  from  the  source  epicenter  of  0 
to  1.95  km  at  an  interval  of  0.05  km,  yielding  a  total  of  40  observation  points.  A  simulated 
fracture  zone  is  located  in  the  third  layer,  and  it  is  shaped  like  a  rectangular  prism  with 
dimensions  of  0.315  km,  0.099  km  and  0.018  km  in  the  i,  y  and  z  directions,  respectively. 
These  axis  directions  are  chosen  such  that  the  receiver  array  is  contained  within  the  i  —  c 
plane  (Figure  8).  We  applied  a  center  frequency  of  25  Hz  for  the  source  wavelet,  and  the 
lattice  spacing  was  set  to  0.009  km  in  all  coordinate  directions.  The  shear  wavelength  is 
9.6  times  the  lattice  spacing  so  that  the  ray-Born  results  should  be  valid  for  back  scattered 
waves.  One  end  of  this  heterogeneity  is  located  under  the  source  point,  and  the  region 
is  centered  under  the  receiver  array  in  the  y  direction,  the  direction  perpendicular  to  the 
receiver  array.  Though  the  background  model  is  one-dimensional,  the  ray-Born  solution 
requires  three-dimensional  calculation  of  the  Green’s  tensors  due  to  the  three-dimensional 
nature  of  the  scattering  lattice.  This  example  allows  the  paraxial  method  to  be  used  to  full 
advantage,  however,  as  rays  need  only  be  traced  down  the  axis  of  the  lattice  in  the  x  —  z 
plane  and  the  paraxial  corrections  can  be  used  to  project  the  results  out  of  this  symmetry 
plane.  Considerable  computation  time  is  saved,  since  the  need  to  trace  rays  along  multiple 
azimuths  is  eliminated. 

To  calculate  the  scattered  field  for  the  isotropically  and  anisotropically  fractured  inhomo¬ 
geneity,  we  use  the  perturbations  to  the  Lame  parameters  given  by  the  Hudson  theory  for  the 
effective  elastic  moduli  of  a  fractured  medium  (Hudson,  1980,  1981;  Crampin,  1984).  The 
effective  moduli  are  given  in  terms  of  an  expansion  to  second-order  in  crack  density  if  =  na^. 
where  n  is  the  number  density  of  cracks,  and  a  is  the  radius  of  the  penny-shaped  cracks. 
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We  applied  a  crack  density  of  0.10,  and  assumed  an  aspect  ratio  of  0.005.  For  the  first  case, 
we  consider  cracks  which  are  dry  or,  equivalently,  filled  with  gas  under  low  pressure.  The 
Lame  parameters  for  the  unfractured  layer  and  the  fractured  perturbations  are  presented  in 
Table  1.  For  the  case  of  perturbations  due  only  to  fine  cracks,  density  variations  are  not 
important  since  the  porosity  volume  is  very  small. 

Shear  Wave  Scattering 

In  the  ray  tracing  for  the  background  field,  we  included  P-  and  S- waves  leaving  the  source 
in  both  the  upward  and  downward  directions.  The  reflections  from  each  interface  were  then 
modeled,  including  phases  with  a  single  mode  conversion  from  P  to  S  or  S  to  P  on  reflection, 
as  waves  changing  mode  more  than  once  are  generally  of  very  low  amplitude.  We  first 
consider  the  effects  of  the  fracture  zone  on  the  pr" ./c-,clion  of  SH  waves  generated  by  a 
single  force  source  oriented  perpendicu'^ir  to  the  receiver  line  (Figure  9).  No  compressional 
waves  are  emitted  in  the  direction  of  the  receiver  array,  and  all  shear  signals  have  a  horizontal 
polarization.  Taking  the  dot  product  of  all  the  Gieen  s  tensors  with  the  source  vector  then 
yields  a  total  of  only  six  arrivals  for  this  source  on  the  transverse  component,  the  shear  to 
shear  wave  reflections  from  each  interface,  including  the  rays  reflecting  first  from  the  free 
surface. 

The  calculated  scattered  field  includes  two  S  to  S  scattered  waves,  one  leaving  the  source 
in  the  downward  direction,  and  the  other  reflecting  from  the  free  surface  and  subsequently 
traveling  to  the  fractured  zone.  The  background  field  and  total  field  (the  sum  of  background 
and  scattered  waves)  for  the  isotropic  fracture  zone  are  compared  in  Figure  10  for  the  time 
interval  noted  in  Figure  9.  Like  the  case  of  reflection  from  a  planar  interface,  the  incident  SH- 
waves  do  not  yield  a  scattered  P-wave  signal  and  only  a  transverse  component  of  displacement 
is  generated.  The  total  displacement  field  distinctly  shows  the  presence  of  the  fractured  zone. 
.Arrows  in  Figure  10  mark  the  arrival  times  of  the  two  scattered  signals  on  the  zero  offset  trace. 


and  they  are  clearly  of  comparable  amplitude  to  the  standard  reflections  from  the  nearby 
planar  interfaces.  The  finite  extent  of  the  fracture  zone  can  be  inferred  from  the  fact  that 
the  back  scattered  waves  are  only  observed  in  the  shorter  offset  receivers.  These  synthetic 
results  show  that  even  thin  fractured  regions  can  be  detected  with  seismic  experiments. 

Next  we  consider  the  case  where  the  experimental  configuration  is  exactly  the  same  as 
above,  except  that  within  the  fracture  zone,  all  the  fractures  are  vertical  and  parallel.  In 
this  anisotropic  case,  the  fractured  zone  is  transversely  isotropic  with  a  horizontal  axis  of 
symmetry.  This  symmetry  direction  provides  an  additional  degree  of  freedom  compared  to 
the  isotropic  case  considered  above,  as  the  orientation  of  the  fractures  relative  to  the  receiver 
array  must  be  considered.  Figures  11  A.  IIB,  and  llC  compare  the  transverse  component 
total  displacement  fields  computed  for  vertical  fractures  aligned  perpendicular,  parallel,  and 
at  45®  to  the  receiver  array,  respectively.  The  resulting  perturbations  to  the  elastic  constants 
are  given  in  Tables  2,  3  and  4.  These  constants  are  different  for  each  case  due  to  the  rotation 
of  the  anisotropic  system  as  the  coordinate  system  is  held  fixed  in  direction.  Note  that  for 
cracks  perpendicular  to  the  receiver  array  6C44  =  0  (Table  2),  which  significantly  affects  the 
total  displacement  field  due  to  the  SH-wave  source.  Since  the  only  incident  strain  component 
which  is  not  vanishingly  small  is  €237  and  the  only  perturbation  in  this  symmetry  system 
which  influences  scattering  in  the  Born  approximation  (equation  7)  is  6C44,  the  scattered 
field  will  be  essentially  zero.  .Accordingly,  in  Figure  11  A,  no  scattered  field  is  present. 

When  the  fractures  are  aligned  parallel  to  the  receiver  array,  the  perturbations  to  the 
clastic  constants  have  the  same  values  as  in  the  previous  anisotropic  model  but  are  rear¬ 
ranged  somewhat  due  to  the  rotation  of  the  elastic  tensor  (Table  3).  The  polarization  of 
incident  SH-waves  is  now  perpendicular  to  the  fracture  plane.  Therefore,  there  will  be  large 
scattered  fields,  as  the  incident  C23  strain  now  interacts  with  a  non-zero  perturbation  (equa¬ 
tion  7).  These  effects  are  seen  in  the  total  field  synthetic  seismograms  for  the  horizontal 
point  force  source  (Figure  1  IB).  Comparison  to  Figure  lOB  shows  that  the  scattered  field 
for  this  anisotropic  model  is  even  larger  than  the  original  results  presented  for  the  isotropic 
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case.  These  results  are  easily  understood  intuitively  as  the  incident  and  scattered  waves 
have  polarization  perpendicular  to  all  of  the  thin  cracks.  Thus  the  maximum  effect  is  to  be 
expected.  This  show.s  the  sensitivity  of  shear  w.avcs  to  fracture  alignment. 

The  results  for  the  models  with  cracks  perpendicular  and  parallel  to  the  receiver  array  are 
special  cases  in  that  both  contain  fractures  oriented  such  that  incident  SH-wave  displacement 
is  polarized  in  a  symmetry  direction  of  the  anisotropic  fracture  zone.  In  the  first  case, 
fractures  perpendicular  to  the  array,  SH  polarization  is  contained  in  the  symmetry  plane. 
In  the  second  example,  the  SH-wave  is  polarized  parallel  to  the  symmetry  axis.  If  the 
fracture  orientation  is  at  some  arbitrary  angle  to  the  receiver  line,  we  would  expect  the 
resulting  synthetic  seismograms  to  be  less  simple.  This  effect  is  shown  by  the  model  where 
the  fractures  are  oriented  at  45°  to  the  array,  with  perturbations  to  the  elastic  constants 
presented  in  Table  4. 

The  transverse  component  total  displacement  field  from  the  horizontal  point  force  source 
is  presented  in  Figure  IIC.  While  this  signal  is  very  similar  to  that  observed  for  the  isotropic 
model  (Figure  lOB)  and  for  the  cracks  parallel  to  the  array  (Figure  IIB),  for  the  first 
time,  the  SH-wave  source  results  in  a  significant  radial  component  synthetic  seismogram 
(Figure  12).  In  the  isotropic  case,  the  perturbations  are  such  that  6C\\  =  8C22  =  8C33  = 

+  Sn,  SC12  =  8Ct3  =  SC23  =  <5^,  and  SC44  =  SC55  =  SCee  =  are  the  only  non¬ 
zero  perturbations.  However,  the  perturbations  for  the  45°  ccise  clearly  have  other  non-zero 
perturbations  (Table  4),  and  in  addition,  the  equivalencies  which  hold  in  the  isotropic  model 
are  not  all  true.  The  result  is  a  complicated  secondary  source  representation  which  yields 
both  SH  and  SV  energy  from  an  incident  SH-wave  (equation  7).  The  vertical  component 
synthetic  seismograms  are  still  essentially  zero  since  the  vertically  propagating  shear  waves 
have  almost  no  vertical  components  of  displacement. 

By  comparing  the  results  for  the  different  orientations  of  the  vertical  cracks  with  respect 
to  the  receiver  array  (Figures  1 1  and  12),  it  is  clear  that  scattered  or  reflected  shear  waves 
from  anisotropic  fracture  zones  are  highly  sensitive  to  the  direction  of  the  displacement 
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polarization  relative  to  the  crack  orientation.  The  returning  energy  from  the  fracture  zone 
varies  in  amplitude  from  a  very  large  signal  to  vanishingly  small,  and  also  shows  changes  in 
polarization  in  intermediate  directions. 

Compressional  Wave  Scattering 

The  utility  of  shear  wave  observations  in  the  inference  of  fracture  zone  anisotropy  is  clear. 
However,  it  is  also  of  interest  to  understand  the  information  given  by  compressional  wave 
experiments  since  it  is  often  easier  to  obtain  good  quality,  high  frequency  P-wave  data  in 
practical  situations.  The  effects  of  the  isotropic  fracture  zone  perturbations  on  P-waves  are 
examined  by  applying  an  explosion  source  to  the  Green’s  tensors.  Figure  13  presents  the 
resulting  background  displacement  field,  which  now  contains  non-zero  radial  and  vertical 
components.  The  largest  disturbances  on  the  vertical  component  are  P-waves,  with  a  pair 
of  reflections  from  the  first  interface  observed  at  0.2  and  0.32  sec  on  the  zero-offset  trace, 
from  the  second  interface  at  0.52  and  0.65  sec,  and  from  the  third  at  0.77  and  0.9  sec.  In 
addition,  some  large  amplitude  shear  arrivals  are  observed  at  larger  offsets.  The  radial  com¬ 
ponent  is  weaker  at  near  source  offsets  since  near  vertical  P-waves  have  minimal  horizontal 
displacement  and  shear  wave  conversions  only  become  significant  at  larger  offsets. 

For  the  gas-filled  fracture  case,  the  perturbation  to  the  Lame  parameter  6X  is  about  5.5 
times  greater  than  the  change  in  rigidity  fi  (Table  1).  VVe  thus  expect  that  the  compressional 
wave  scattering  will  be  relatively  significant,  as  a  change  in  A  affects  only  P-wave  scattering 
(Gibson  and  Ben-Menahem,  1991).  Since  the  nearly  vertically  travelling  compressional  waves 
have  a  dominantly  vertical  displacement,  only  the  vertical  component  synthetic  seismograms 
for  the  total  displacement  field  for  the  explosion  source  are  shown  in  Figure  14.  where  the 
arrival  times  of  the  principal  scattered  waves  are  marked.  The  scattered  compressional 
waves  are  easily  observed  on  the  vertical  component.  These  results  and  the  example  of  the 
transverse  component  SH-waves  show  that  both  compressional  and  shear  waves  are  of  use 
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in  the  detection  of  isotropic,  gas-filled  fracture  zones. 

In  contrast,  the  scattering  of  incident  compressional  wave  energy  by  anisotropic  fracture 
zones  yields  almost  no  information  on  crack  orientation.  Figure  15  displays  the  total  field 
synthetic  seismograms  for  the  three  models  considered  for  the  shear  waves,  with  fractures 
aligned  perpendicular,  parallel,  and  at  45°  to  the  receiver  array.  Comparison  of  the  results 
for  the  three  models  shows  that  the  predicted  synthetic  seismograms  are  essentially  the 
same,  and  that  the  properties  of  the  scattered  P-waves  do  not  depend  on  the  azimuthal 
orientation  of  vertical  fracture  zones.  Mathematically,  this  can  be  easily  understood  as  the 
result  of  the  strains  induced  by  the  incident  wavcficld  and  the  consequent  secondary  source 
terms  due  to  elastic  constant  perturbations  appearing  in  equation  (7).  The  incident  P- 
wave  strain,  mostly  £33,  interact  principally  with  the  perturbation  6C33.  Though  reduced 
from  the  isotropic  equivalent  <5A  -I-  28 (Table  1),  this  perturbation  is  still  large  enough  to 
generate  energy  on  the  vertical  component  seismograms  so  that  the  overall  trends  displayed 
by  the  scattered  waves  are  similar  to  those  for  the  isotropic  fracture  zone  (Figure  14B), 
though  slightly  weaker.  However,  rotation  of  the  vertical  fractures  with  respect  to  the 
receiver  array  leaves  ^033  invariant  (compare  Tables  2,  3  and  4).  Therefore,  the  scattered 
waves  are  practically  constant  with  fracture  orientation  and  we  see  that  at  least  for  nearly 
vertically  incident  P-waves,  110  information  on  fracture  orientation  is  available.  The  physical 
interpretation  of  these  predictions  is  also  straightforward,  as  the  displacement  of  the  incident 
wavefield  is  vertical  and  therefore  more  or  less  parallel  to  the  crack  plane  for  all  orientations 
of  cracks:  we  should  not  expect  much  variation  of  reflection  properties.  This  also  explains  the 
reduced  amplitude  of  the  scattered  waves  from  the  anisotropic  zone  (Figure  15)  as  compared 
to  the  isotropic  fractures  (Figure  14). 
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APPLICATION  TO  FIELD  DATA 


The  delineation  of  fracture  zones  is  of  great  practical  interest  in  the  development  of  geother- 
maJ  fields,  as  the  permeability  created  by  the  fractures  can  control  the  fluid  flow  in  the 
geothermal  systems.  To  facilitate  the  development  of  the  Larderello  geothermal  field  in  ItaJy 
(Figure  16),  an  extensive  geophysical  study  has  been  conducted  (Batini  et  ai,  1983,  1985a. 
1985b,  1985c,  1991).  Recently,  a  VSP  experiment  was  conducted  in  the  Badia  lA  well  in 
the  Larderello  field  to  attempt  to  further  delineate  several  possible  fracture  zones  (Batini  et 
ai,  1991).  Interpretation  of  the  processed  data  from  the  VSP  experiments  and  a  number 
of  other  surface  seismic  surveys  resulted  in  the  suggestion  that  three  reflected  events  from  a 
subsurface  feature  below  the  depth  of  the  well  bottom  and  labeled  the  “H  marker”  were  due 
to  fracture  zones  distributed  over  a  depth  range  of  approximately  0.200  km.  This  model  is 
based  on  some  regional  seismic  observations.  We  have  used  ray- Born  synthetic  seismograms 
to  test  this  a  priori  hypothesis  and  to  develop  a  laterally  varying  model  of  fracture  density 
in  three  fracture  zones,  which  accounts  for  these  arrivals  as  observed  from  two  VSP  surveys. 

Background  Model 

A  local  three-dimensional  seismic  survey  allowed  the  development  of  a  fairly  detailed  model 
of  the  principal  geological  layers  in  the  Badia  area  (Batini  et  ai,  1991).  The  principal 
geological  features  of  the  locality  are  two  shallow  layers  consisting  of  various  sedimentary 
units  overlying  a  thicker  zone  of  metamorphic  rock.  The  LI  and  L  markers  form  the  lower 
interfac<-s  of  the  first  and  .second  sedimentary  layers,  respectively.  P-wave  velocities  for  the 
model  are  given  in  a  cross-section  in  Figure  17.  a  section  in  the  east- west  direction  which 
intersects  the  position  of  the  Badia  lA  well.  Contour  maps  of  the  depths  of  the  LI,  L,  and 
H  markers  arc  presented  in  Figure  18,  along  with  the  positions  of  the  Badia  lA  well  and  the 
two  locations  for  the  Vibroseis  source.  One  of  tJie  shot  points  is  very  close  to  the  well  and 
will  be  called  the  zero  offset  point,  and  the  A  shot  point  is  0.981  km  north  and  slightly  east 


of  the  well.  These  maps  clearly  demonstrate  the  three-dimensional  nature  of  this  modeling 
problem. 

Data 

Due  to  difficult  drilling  conditions,  the  Badia  lA  well  is  highly  deviated  (Figure  19),  so  that 
only  vertical  gcophone  components  have  sufficiently  high  signal-to-noise  ratio  to  analyze  the 
comparatively  weak  scattered  arrivals.  Many  of  the  vertical  component  data  traces  were, 
however,  extremely  noisy  after  the  first  arrival,  especially  at  shallower  depths.  Therefore, 
a  subset  of  recorded  data  was  used,  with  receiver  positions  indicated  in  Figure  19.  Slightly 
different  depth  ranges  are  covered  by  the  geophones  from  the  zero  and  A  offset  experiments. 

Data  from  both  shot  points  are  shown  in  Figure  20  where  each  trace  is  normalized  to 
unit  amplitude  on  the  first  arrival.  The  signal  from  the  far  offset  VSP  is  somewhat  noisier, 
since  the  first  arrival  prior  to  normalization  is  lower  in  amplitude  due  to  a  greater  distance 
from  source  to  receiver.  Therefore,  noise  is  amplified  more  in  the  normalization  than  with 
the  zero-offset  traces.  For  both  data  sets,  however,  the  downgoing  waves  are  clearly  much 
stronger  than  any  reflected  arrivals,  which  are  not  visible  in  these  sections.  To  bring  out 
the  upgoing  reflected  waves,  we  processed  the  data  using  a  median  filter  (Hardage,  1983; 
Reiter.  1991).  The  moveout  velocity  on  the  filter  was  set  to  the  opposite  of  the  apparent 
velocity  of  the  downgoing  wave  to  enhance  the  reflected  P-waves,  and  the  filter  was  applied 
across  17  traces.  After  median  filtering,  the  data  were  subsequently  low-pass  filtered  with 
cutoff  frequencies  of  5.5  and  85  Hz  for  the  zero  and  A  offset  data,  respectively.  Since  the 
predominant  signal  strength  is  at  30  Hz,  this  should  not  significantly  affect  the  desired  signal 
in  the  results.  The  processed  data  are  displayed  in  Figure  21  with  a  time  squared  gain  factor 
applied,  and  the  zero  offset  data  are  magnified  by  a  factor  of  three  compared  to  the  A  offset 
plot.  The  arrivals  from  the  H  marker  are  indicated  on  these  plots.  Considering  the  difference 
in  gain  factors  applied  to  the  upgoing  data  plots,  it  is  clear  that  the  signal  observed  from 
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the  A  offset  source  is  much  stronger,  suggesting  a  lateral  difference  in  the  properties  of  the 
II  marker. 


Modeling  Results 

The  geological  model  (Figures  17  and  18)  clearly  demonstrates  the  three-dimensional  struc¬ 
ture  of  the  sedimentary  and  metamorphic  layers  in  the  Larderello  area.  This  significantly 
affects  the  ray  tracing  procedures  which  must  illuminate  the  H  marker  to  calculate  the  scat¬ 
tered  waves,  as  would  be  expected.  To  illustrate  this  problem,  we  show  the  ray  paths  for 
a  fan  of  rays  traced  from  both  the  zero  and  A  offset  source  positions  with  azimuthal  take¬ 
off  angle  in  the  north-south  direction  (Figure  22).  The  resulting  ray  paths  for  both  source 
positions  have  large  distortion  in  the  east- west  direction  due  to  the  Ll  and  L  interfaces  (Fig¬ 
ure  18),  whereas  in  a  one-dimensional  earth  model,  these  rays  would  be  entirely  contained  in 
a  vertical  plane.  In  particular,  we  see  that  for  the  A  offset  source  the  bending  of  the  rays  at 
these  interfaces  is  so  large  that  the  rays  exit  the  model  before  reaching  the  H  marker  depth. 
Therefore,  computation  of  the  Green’s  tensors  for  the  incident  wavefield  in  equation  (6)  re¬ 
quires  that  we  trace  fans  of  rays  like  those  in  Figure  22  over  all  azimuths.  Approximately 
2,800  rays  were  traced  from  the  two  source  positions.  Although  the  ray  tracing  is  simpler 
from  the  receivers  to  the  H  marker  since  the  metamorphic  layer  is  homogeneous,  even  more 
rays  are  required  since  the  closer  proximity  of  the  receivers  to  the  marker  increases  the  range 
of  take-off  angles  for  the  fans  of  rays  over  all  azimuths.  Almost  8,000  rays  were  traced 
from  each  receiver  to  insure  sufficient  coverage  of  the  H  marker.  These  points  regarding 
the  tlircxi-dimensional  nature  of  the  wave  propagation  for  this  problem  emphasize  the  value 
of  the  ray  tracing  approach,  which  is  a  practical  and  feasible  method  of  solution.  A  more 
complete  solution  might  in  principle  be  obtained  using  a  finite-difference  method,  but  the 
required  discretization  of  the  three-dimensional  model  to  adequately  represent  the  different 
layers  and  the  H  marker  would  be  prohibitively  time  consuming. 
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Since  the  scattered  waveforms  are  stronger  on  the  A  offset  data,  we  began  by  developing 
a  model  which  would  account  for  these  data.  A  lattice  was  set  up  to  conform  to  the  map 
of  the  H  horizon  in  Figure  18C,  though  due  to  limitations  of  computer  storage,  we  limited 
the  size  of  the  lattice  to  1.20  by  1.20  km.  It  is  more  efficient  for  this  modeling  of  scattering 
by  a  thin  sheet  to  actually  restrict  the  lattice  points  to  trace  the  depth  of  the  sheet  rather 
than  specifying  a  regular  three-dimensional  Cartesian  lattice  which  would  have  many  zero 
valued  nodes.  For  the  A  shot  point,  the  lattice  ranged  from  0.6  to  1.8  km  in  the  east- west 
direction,  the  complete  width  of  the  map,  but  only  from  0.5  to  1.7  km  in  the  north-south 
direction.  Similarly,  the  lattice  for  the  zero  offset  source  ranged  from  O.S  to  2.0  in  the  north- 
south  direction.  Utilizing  a  discretization  in  the  two  horizontal  directions  of  0.020  km.  this 
results  in  a  three-dimensional  lattice  of  61  by  61  by  1  points  mapping  the  H  marker.  Other 
forms  of  data  analysis  were  unable  to  specify  a  thickness  for  the  hypothetical  fracture  zones, 
though  various  estimates  ranged  from  0.010  to  0.060  km.  We  therefore  set  the  thickness  of 
the  lattice  to  be  0.020  km.  A  source  pulse  for  both  the  zero  and  A  offsets  wcis  applied  by 
choosing  the  waveform  recorded  as  a  first  arrival  on  representative  traces  with  good  signal- 
to-noise  ratio.  Final  synthetic  seismograms  should  therefore  be  directly  comparable  to  the 
field  observations. 

By  trial  and  error  forward  modeling,  it  was  determined  that  a  model  consisting  of  three 
fractured  horizons  could  account  for  the  observations.  This  was  accomplished  by  temporarily 
neglecting  amplitude  effects  and  matching  only  the  arrival  times  of  the  observed  data.  The 
depth  and  shape  of  the  first  of  the  zones,  the  HI  event,  was  left  to  conform  to  the  map  in 
Figure  18C,  since  it  was  obtained  by  the  regional  three-dimensional  survey.  The  second  and 
third  zones,  the  H2  and  H3  events,  were  defined  to  have  the  same  lateral  variation  as  shown 
in  Figure  ISC,  but  deeper  by  0.090  and  0.170  km,  respectively.  Our  analysis,  therefore, 
confirms  a  total  thickness  of  the  H  marker  on  the  order  of  0.200  km.  It  should  be  noted  that 
the  largest  amplitudes  of  the  scattered  waves  in  Figure  21  correspond  to  the  arrivals  from 
the  H2  and  H3  horizons,  in  agreement  with  some  other  regional  observations  (Batini  et  al.. 
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1991). 

After  using  these  initial  stages  of  modeling  to  define  the  depths  of  the  HI.  H2  and  H3 
markers,  we  developed  a  model  to  account  for  the  differing  strengths  of  the  recorded  scattered 
waves  from  the  two  offsets  by  considering  the  regions  imaged  by  the  two  experiments.  In 
Figure  23,  we  present  contour  maps  of  the  tot2d  travel  time  from  source  to  lattice  point  on 
the  H2  marker  and  back  to  the  receiver  for  two  pairs  of  receivers  at  equivalent  well  positions. 
The  dot  in  the  interior  of  each  plot  indicates  the  minimum  time,  which  corresponds  to  the 
Fermat’s  principle  true  travel  time  for  the  reflected  pulse.  The  contour  interval  for  these  plots 
was  chosen  so  that  each  contour  outlines  a  Fresnel  zone.  Sheriff  and  Geldart  ( 19S2)  define  the 
Fresnel  zone  such  that  the  distance  from  source  to  scattering  surface  is  a  quarter  wavelength 
larger  for  the  outer  edge  of  the  zone  than  from  the  inner  edge.  Using  the  dominant  frequency 
for  these  data  of  30  Hz  and  the  medium  velocity  of  5.1  km/sec,  this  distance  corresponds  to  a 
travel  time  difference  of  0.0083  sec.  Since  most  of  the  effect  of  the  scattering  surface  is  from 
that  region  enclosed  within  the  first  Fresnel  zone,  these  figures  allow  a  simple  determination 
of  the  region  of  the  H2  marker  observed  in  the  two  VSP  experiments. 

Comparison  of  Figures  23A  and  23C  for  receivers  at  the  shallow  region  of  the  data  and 
Figures  23B  and  23D  from  the  deep  portion  shows  that  the  A  offset  survey  is  imaging  a 
region  of  the  H  marker  several  hundred  meters  further  to  the  north  than  the  zero  offset  shot 
point,  though  the  imaging  point  for  the  deeper  A  offset  receiver  (Figure  23D)  approaches 
that  for  the  shallow  zero  offset  receiver  (Figure  23A).  We  developed  a  model  assuming 
that  the  velocity  variations  at  these  markers  are  due  to  fracturing  and  that  there  must  be 
some  larger  velocity  perturbations  toward  the  north,  toward  shot  point  A.  .^fter  attempting 
several  models,  an  isotropic  model  was  developed  which  satisfactorily  explains  the  data.  The 
crack  density  in  the  HI  marker,  which  is  weak  in  both  data  sets  was  set  to  0.004  uniformly, 
whereas  for  the  H2  and  113  markers,  crack  density  was  set  to  0.14  north  of  the  1.4  km  line, 
and  0.004  south  of  this  border.  Temperature  and  pressure  conditions  at  the  depth  of  the  H 
marker  are  very  uncertain,  but  other  wells  suggest  that  values  of  about  300°C  and  several 
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hundred  bar  pressure  are  appropriate  (Batini  et  ai,  1983).  Under  these  conditions,  the  bulk 
modulus  of  the  pore  filling  fluid  can  probably  be  roughly  estimated  as  0.1  GPa  (Anderson 
and  Whitcomb,  1973).  Table  6  shows  the  resulting  perturbations  to  Lame  parameters.  The 
synthetic  seismograms,  processed  e.xactly  as  were  the  data,  are  shown  in  Figure  24.  From 
these  results  we  infer  that  our  model  adequately  accounts  for  most  of  the  seismic  properties 
in  the  region  imaged  by  the  A  source.  The  quality  of  the  fit  of  the  zero  offset  synthedetic  is 
less  satisfactory,  however,  which  is  at  least  in  part  a  consequence  of  the  lower  signal  to  noise 
ratio  for  the  weak  scattered  wave  in  this  data  set.  At  least  part  of  the  signal,  especially  from 
traces  5-5  to  65,  is  still  relatively  well  modeled. 

Although  some  observations  from  other  wells  suggested  that  vertical  fractures  are  present 
in  the  H  marker,  the  successfi  l  match  of  data  and  synthetic  seismograms,  particularly  for 
the  A  offset  data,  suggest  that  a  purely  isotropic  model  can  account  for  these  seismic  obser¬ 
vations.  The  synthetic  examples  shown  above  for  the  explosion  source  clearly  demonstrate 
that  in  any  r,ase  a  compressional  wave  experiment  will  not  yield  any  information  allowing 
a  unique  interpretation  of  the  presence  of  vertical  fractures,  the  only  difference  from  the 
isotropic  case  being  a  weaker  reflected  wave.  We  attempted  to  develop  a  model  for  the  H 
marker  using  the  eltistic  constant  perturbations  appropriate  for  vertical  fractures,  but  even 
a  crack  density  as  high  as  0.30  yielded  a  signal  for  the  A  offset  synthetics,  which  was  far 
weaker  than  the  observed  data.  This  crack  density  is  already  so  high  that  it  likely  violates 
the  single  scattering  assumption  used  in  the  derivation  of  the  expressions  for  the  effective 
elastic  constants,  so  we  did  not  test  any  higher  values.  Instead,  we  interpret  these  results 
as  suggesting  that  the  principal  source  of  the  scattered  seismic  energy  is  an  isotropic  veloc¬ 
ity  change  which  may  have  some  weak  anisotropy  superimposed  due  to  vertical  fractures. 
The  available  data  does  not  allow  any  further  conclusions,  and.  as  shown  in  the  synthetic 
examples,  shear  wave  observations  would  be  necessary  to  uniquely  determine  the  presence 
of  anisotropy  due  to  vertical  fractures. 
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DISCUSSION  AND  CONCLUSIONS 


The  ray-Born  method  for  the  modeling  of  the  effects  of  small,  localized  regions  of  inhomo¬ 
geneity  on  seismograms  is  general  and  can  in  principle  be  applied  to  fully  anisotropic  models. 
Limitations  on  the  ray  theoretical  aspects  of  the  algorithm  result  from  the  well  known  as¬ 
sumption  of  high  frequency  methods,  namely,  that  the  wavelength  of  the  propagating  signal 
must  be  much  less  than  the  scale  length  of  the  background  earth  model  (Ben-Menahem  and 
Beydoun,  1985).  Comparisons  of  ray-Born  solutions  with  the  complete,  discrete  wavenumber 
solution  for  the  elastic  wave  scattering  from  a  spherical  inhomogeneity  illustrate  some  of  the 
limitations  of  the  Born  approximation.  The  principal  restriction  is  also  one  of  scale  lengths. 
As  the  ratio  of  the  propagating  wavelength  to  the  scale  length  of  the  inhomogeneity  decreases 
below  a  value  of  about  5,  the  Born  approximation  fails  to  reproduce  some  significant  features 
of  the  scattered  waves.  For  the  sphere,  the  missing  features  were  the  dual  reflections  from 
the  front  and  back  of  the  sphere  in  the  back  scattering  direction  (Figures  5  and  6).  The 
ray-Born  solution  predicted  only  a  single,  broad  pulse. 

With  these  limitations  in  mind,  we  applied  the  method  to  synthetic  studies  of  a  thin, 
laterally  extended  fracture  zone  in  a  simple  layered  earth  model  (Figure  8).  These  synthetic 
results  clearly  demonstrate  that  we  expect  shear  waves  to  provide  useful  information  on 
the  alignment  direction  of  vertical  fractures,  as  the  scattered  wavefield  varies  significantly 
with  the  polarization  direction  of  the  incident  shear  wave.  In  particular,  if  the  incident  SH- 
wave  is  polarized  at  an  angle  to  the  fracture  orientation,  the  scattered  wavefield  can  have 
significant  energy  on  the  radial  component.  Similar  conclusions  were  reached  by  Spencer  and 
Chi  (1991)  in  a  theoretical  examination  of  vertically  incident  shear  waves  on  a  uniformly 
fractured  layer  or  half-space.  Mueller  (1991)  showed  the  correlation  of  lateral  variation 
of  reflection  properties  of  SH  and  SV  waves  with  fracture  intensity  in  the  Austin  Ch.alk 
located  in  (,'entral  Texas.  These  effects,  especially  when  combined  with  the  variation  of  the 
observations  witli  different  incident  polarizations,  are  not  likely  to  be  reproduced  by  realistic 
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isotropic  iiihomogcncity  in  similar  earth  models. 

The  synthetic  models  also  show  that  compressional  wave  scattering  varies  little  with 
crack  orientation,  although  we  can  still  use  reflected  waves  to  detect  fractured  regions.  A 
more  unique  aspect  of  the  P-wave  data  is  its  sensitivity  to  the  material  filling  the  cracks. 
Since  perturbations  to  the  isotropic  Lame  parameter  A  affect  only  P-wave  scattering  and 
the  presence  of  liquid  in  the  cracks  decreases  SX  significantly  from  the  gas-filled  case,  the 
amplitude  of  reflected  P- waves  also  decreases  dramatically  (Gibson.  1991).  These  ray-Born 
results  reproduce  the  classic  ‘‘bright  spot’’  behavior  often  used  to  distinguish  gas  reservoirs 
from  oil  reservoirs. 

These  guidelines  help  to  interpret  the  VSP  data  from  the  Badia  lA  well  in  the  Italian 
Larderello  geothermal  field.  VVe  developed  an  isotropic  model  of  three  fracture  zones  con¬ 
tributing  to  the  seismic  H  marker,  accounting  for  observed  reflected  waves  from  below  the 
depth  of  the  well.  The  model  implies  that  fracture  density  increases  northwards  and  that  this 
would  be  a  good  direction  for  further  exploratory  drilling  in  development  of  the  geothermal 
resources. 

There  are,  however,  several  ambiguities  in  the  modeling  which  cannot  be  removed  due 
to  limitations  on  the  available  data  and  which  will  affect  all  seismic  studies  of  fractures, 
especially  those  using  P-wave  data.  For  example,  we  cannot  rule  out  the  presence  of  aligned, 
vertical  fractures  in  the  H  marker  since  P-wave  scattering  from  a  vertically  fractured  region 
shows  no  indication  of  fracture  orientation.  Comparison  of  synthetic  models  for  isotropic 
and  anisotropic  fracture  zones  of  equal  crack  density  (e.g.,  Figures  14  and  15)  shows  that 
the  amplitude  of  the  observed  P-waves  should  be  reduced  for  the  anisotropic  model.  The 
large  scattered  field  in  the  A  offset  data  was  not  reproducible  with  realistic  anisotropic 
crack  densities,  which  gives  an  indirect  evidence  for  non- vertical  fractures  in  this  region.  An 
ambiguity  in  the  ray-Born  modeling  of  fracture  zones  that  enters  here  is  the  state  of  the  pore 
fluids.  VVe  applied  a  value  for  bulk  modulus  of  the  fluid  ba.sed  on  some  values  representative 
for  the  estimated  in  situ  conditions.  It  is  possible,  however,  that  there  is  some  error  in  the 
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hulk  modulus  value.  If  it  is  significantly  overestimated,  the  effects  of  the  fractures  will  be 
more  like  those  of  the  synthetic  gas-filled  fracture  models  with  non-zero  bulk  modulus  of 
pore  fluid  and  P-wave  scattering  amplitudes  will  increase  for  a  given  'rack  density.  This 
might  make  it  possible  to  develop  an  anisotropic  model  accounting  for  ..re  field  data. 

Another  area  of  ambiguity  is  the  equivalent  effects  of  increased  fracture  zone  thickness 
and  increased  fracture  density.  Due  to  the  linearity  of  the  method  and  our  representation 
of  the  scattering  lattices  as  only  a  single  unit  cell  in  the  vertical  direction,  a  doubling  of  the 
thickness  of  the  unit  cell  results  in  a  doubling  of  the  scattered  wavefield  (equation  8).  The 
same  effect  results  from  a  doubling  of  the  perturbations  to  elastic  constants.  Our  assumed 
value  of  0.020  km  for  the  thickness  of  the  three  fractured  intervals  Hi,  II2  and  H3  is  of  the 
same  order  of  size  as  other  estimates  (Batini  et  ah,  1991),  but  could  very  well  be  off.  At 
the  same  time,  the  successful  modeling  of  the  observed  data,  especially  for  source  offset  A. 
suggests  that  our  model  of  the  zones  as  a  single  lattice  point  in  thickness  is  not  far  off  and 
that  the  zones  are  not  too  thick. 

It  is  possible  that  both  isotropic  heterogeneity  and  superimposed  vertical  fracturing  are 
present.  In  this  active  geothermal  area,  ongoing  hydrothermal  processes  will  likely  cause 
mineral  alteration  and  the  sealing  of  many  fractures  (Batini  et  ah.  1983:  Batini  et  ah, 
1985c).  Under  these  conditions,  the  properties  of  the  rock  across  the  fractured  regions 
could  be  altered  in  such  a  way  that  a  superposition  of  anisotropic  fracturing  and  isotropic 
velocity  variations  is  not  unrealistic.  An  increase  in  porosity  due  to  pores  of  large  aspect 
ratio  would  change  density  proportionately  more  the  elastic  properties  and  would  lead  to 
isotropic  scattering  of  elastic  waves.  It  appears  from  the  modeling  of  the  P-wave  data  that 
the  velocity  variations  are  not  entirely  due  to  vertical  fracturing  in  any  case,  and  an  isotropic 
model  can  explain  the  data.  The  only  way  to  concretely  determine  the  presence  or  absence 
of  anisotropy  would  be  to  obtain  high  quality  shear  wave  data  from  the  same  locations. 

Similar  conclusions  regarding  the  information  contained  in  seismic  data  were  obtained  by 
Slolt  nnd  VVeglein  (1985)  in  an  analysis  of  mult iparamcler,  lim‘arize«l  inversion  iiK'thods.  Due 
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to  the  physics  of  wave  propagation  and  limitations  on  observed  quantities  from  experimental 
configurations,  signal-to-noise  ratio,  and  the  deep  exploration  targets,  it  may  only  be  possible 
to  achieve  some  knowledge  of  the  location  in  depth  or  time  of  heterogeneity.  Increasingly 
ambitious  goals  under  better  conditions  or  more  elaborate  experiments  include:  (1)  the 
sign  of  the  property  changes  in  the  heterogeneous  zone;  (2)  the  magnitude  of  the  changes; 
(3)  lateral  variations  of  data  amplitudes  allowdng  more  detailed  analysis;  (4)  inference  of 
multiple  physical  properties;  and  (5)  true  values  of  all  earth  properties.  Clearly  the  last 
goal  is  only  going  to  be  achievable  in  extremely  rare  cases.  The  analysis  of  the  properties 
of  scattered  w’avefields  from  fractured  zones  shows  that  to  achieve  anything  but  the  first 
two  or  three  goals  requires  very  good  quality  seismic  data  and  must  incorporate  shear  wave 
observations.  Therefore,  given  that  only  compressional  wave  data  of  relatively  low  signal  to 
noise  ratio  was  available,  the  ray- Born  method  allowed  some  very  useful  information  on  the 
positions  both  vertically  and  laterally  of  some  changes  in  earth  properties. 

Our  models  of  both  synthetic  and  field  data  emphasize  the  utility  of  the  ray-Born  method 
for  modeling  scattered  wavefields  in  complicated  three-dimensional  geological  structures.  Al¬ 
though  there  is  non-uniqueness  in  relating  amplitude  to  fracture  density  and  fracture  zone 
thickness,  the  model  for  the  Badia  lA  data  does  well  in  predicting  the  kinematic  properties 
of  the  scattered  waves.  It  should  therefore  give  concrete  and  valuable  information  on  the 
depths  of  the  fracture  zones.  This  type  of  information  is  of  great  utility  in  geothermal  field 
development.  Especially  for  large  three-dimensional  problems,  the  method  is  a  compara¬ 
tively  rapid  and  efficient  means  of  exploring  the  effects  of  different  models  of  small  scale 
heterogeneity  within  the  earth.  The  results  also  clearly  show  the  value  of  shear  wave  obser- 
vat  ons  in  detailed  resolution  of  subsurface  seismic  properties,  particularly  for  the  resolution 
of  anisotropy. 
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Parameter 

Background  Value  (GPa) 

Perturbation  (GPa) 

A 

21.49 

-8.34 

11.86 

-1.49 

Table  1.  Background  Lame  parameters  and  perturbations  for  the  fracture  zone  in  the  model 
in  Fig.  8  with  randomly  oriented,  gas-filled  fractvircs. 


Parameter 

Background  Value  (GPa) 

Perturbation  (GPa) 

Cii 

45.21 

-18.6 

C22  =  G33 

45.21 

-4.19 

C23 

21.49 

-4.19 

C'la  =  C'u 

21.49 

-8.80 

C44 

11.86 

0.0 

C55  =  Cae 

11.86 

-2.29 

Table  2.  Background  elastic  constants  and  perturbations  lor  the  fracture  zone  in  the  model 
in  Fig.  8  with  gas-filled  vertical  fractures  oriented  perpendicular  to  the  receiver  array.  The 
parameters  C\\,Ci2i  and  C33  are  all  equivalent  to  A  -|-  Ifi  in  an  isotropic  medium,  while 
C12  =  Ci3  =  C23  =  A  and  C44  =  C55  =  Cee  =  /*• 
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Parameter 

Background  Value  (GPa) 

Perturbation  (GPa) 

Cii  =  C33 

45.21 

-4.19 

C22 

45.21 

-18.6 

Cia 

21.49 

-4.19 

^“23  =  f^l2 

21.49 

-8.80 

C55 

11.86 

0.0 

C44  =  Cr,6 

11.86 

-2.29 

Table  3.  Background  elastic  constants  and  perturbations  for  the  fracture  zone  in  the  model 
in  Fig.  8  with  gas-filled  vertical  fractures  oriented  parallel  to  the  receiver  array. 


Parameter 

Background  Value  (GPa) 

Perturbation  (GPa) 

Cii  =  C22 

45.21 

-12.4 

C33 

45.21 

-4.19 

Gis  =  C23 

21.49 

-6.50 

C12 

21.49 

-7.S0 

C44  =  (I'ss 

11.86 

-1.15 

Gee 

11.86 

-1.29 

Gi6  =  G26 

0.0 

3.60 

G36 

0.0 

2.31 

C4, 

0.0 

1.15 

Table  4.  Background  elastic  constants  and  perturbations  for  the  fracture  zone  in  the  model  in 
Fig.  8  with  gas-filled  vertical  fractures  oriented  at  15®  to  the  receiver  array.  In  an  isotropic 
medium  or  in  a  transversely  isotropic  material  in  a  coordinate  system  where  one  of  the 
coordinate  axes  parallels  the  axis  of  symmetry,  Cio,  C261C  36-  and  C'^s  are  all  zero. 
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Crack  density 

Parameter 

Background  Value  (GPa) 

Perturbation  (GPa) 

^  =  0.14 

A 

23.6 

-3.33 

23.6 

-3.33 

/>* 

II 

o 

o 

o 

A 

23.6 

-0.10 

1 

23.6 

-0.10 

Table  5.  Background  Lame  parameters  and  perturbations  for  the  fracture  zone  model  for 
the  Badia  location  H  markers.  The  marker  HI  marker  was  assigned  a  crack  density  of  0.004 
throughout  the  model,  whereas  the  computed  synthetics  applied  a  crack  density  of  0.14  to 
the  H2  and  H3  markers  north  of  the  1.4  km  latitudinal  line  (Fig.  18). 
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Figure  1:  Discretization  scheme  for  implementation  of  the  ray- Born  algorithm.  A  two- 
dimensional  example  is  shown  for  clarity,  though  the  actual  implementation  includes 
three-dimensional  scattering  lattices  and  background  models.  All  lattice  points  outside 
the  boundary  of  the  inhomogeneity  are  assigned  zero  perturbation  values,  while  those 
points  inside  the  boundary  can  be  assigned  non-zero  perturbations  to  density  or  any 
of  the  elastic  constants  C/j.  As  the  incident  wavefront  encounters  a  lattice  point,  that 
particular  point  becomes  a  secondary  source. 
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Figure  2:  Diagrammatic  illustration  of  the  ray  coordinates  71  and  72  and  ray-centered  co¬ 
ordinate  basis  vectors  6],  62,  and  63.  The  ray  coordinates  correspond  to  the  azimuthal 
and  declination  take-off  angles,  71  and  72,  respectively.  At  each  point  P  on  the  ray,  the 
two  basis  vectors  ei  and  62  are  located  in  the  plane  tangent  to  the  wavefront,  while  the 
third  vector  63  is  tangent  to  the  ray,  forming  a  right-handed  coordinate  system  valid  in 
the  vicinity  01  the  ray  (Cerveny,  1985). 
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Figure  3:  Sphericad  coordinate  f.ystem.  The  scattered  wavefields  from  the  elastic  sphere  are 
calculated  for  receiver  distances  sufficiently  distant  that  the  P-wave  displacement  will  be 
observed  entirely  on  the  radial  component  e,  and  the  S-wave  will  appear  only  on  the 
eo  component.  The  shear  wave  has  only  a  ej  component  due  to  the  symmetry  of  the 
problem  for  a  plane  P-wave  incident  along  the  z  axis. 


46 


-60  -40  -20  0  20  40  60 

X  -COORDINATE  (KM) 


Figure  4:  Receiver  configuration  for  the  computation  of  wavefields  scattered  from  an  elastic 
sphere.  The  sphere  is  located  at  the  origin  of  the  coordinate  system  and  the  incident 
plane  P-wave  is  propagating  in  the  downwards  vertical  direction  so  that  the  direction 
90°  corresponds  to  forward  scattering,  and  270°  is  the  direction  for  back  scattered  waves. 
The  angles  given  here  are  equal  to  the  angular  coordinate  6  in  Figure  3  minus  90°  and 
are  used  for  reference  in  the  synthetic  seismogram  plots. 
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Figures  5A  and  5B:  ray-Born  radial  component  synthetic  seismograms  for  the  scattering 
from  a  compressional  plane  wave  vertically  incident  on  a  sphere  showing  scattered  com- 
pressional  waves.  Since  the  amplitudes  of  the  scattered  waves  increase  dramatically  with 
increasing  frequency,  the  scale  of  the  different  plots  was  changed  to  allow  the  different 
waveforms  to  be  observed.  A)  Ratio  of  compressional  wavelength  to  sphere  diameter 
ijp  =  18,  frequency  /  =  0.25  Hz.  B)  7p  =  I..5,  /  =  1  Hz.  Blotting  scale  multiplied  by  0.1 
relative  to  the  rjp  =  18  plot. 
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Figures  oC  and  5D:  C):  rjp  =  l.S,/  =  2.5  Hz.  Plotting  scale  multiplied  by  0.02  relative  to 
the  7/p  =  IS  plot.  D)  7/p  =  0.!),/  =  5  Hz.  Plotting  scale  multiplied  by  0.01  relative  to  the 
7/p  =  IS  plot. 
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Figures  6A  and  6B:  Discrete  wavenumber  radial  component  synthetic  seismograms  for  the 
scattering  from  a  compressional  plane  wave  vertically  incident  on  a  sphere  showing  scat¬ 
tered  compressional  waves.  All  plot  scales  for  a  given  frequency  are  the  same  as  in  Figure  5 
so  that  waveforms  may  be  directly  compared  for  the  ray-liorn  and  discrete  wavenumber 
methods.  A)  Ratio  of  compressional  wavelength  to  s[«here  diameter  rjp  =  18.  frequency 
/  =  0.25  Hz.  B)  rip  =  4.5,  /  =  1  Hz. 
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Figures  60  and  6D:  C):  7/p  =  1.8.  /  =  2.5  Hz.  D)  7/p  =  0.9,/  =  5  Hz. 
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I'  igun;  7;  Comparison  of  the  arnpliliules  of  coinpressional  waves  scattered  from  tlie  sphere 
with  1%  velocity  perturbations  with  different  values  of  incident  wavelength  to  sphere 
diameter  ratio  r/p  indicated.  Forward  scattering  is  in  the  direction  90"  and  back  scattering 
is  at  270". 
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Figure  8:  Background  model  used  to  calculate  synthetic  seismograms  including  the  effect 
of  a  thin  fracture  zone.  Receivers  are  indicated  by  the  black  diamond  symbols,  and  the 
position  of  the  source  is  shoAvn.  This  model  cross-section  is  contained  in  the  x  —  z  plane, 
and  the  y  axis  is  perpendicular  to  the  figure. 
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Figure  9;  Transverse  component  background  displacement  field  from  a  cross-line  point  source 
for  the  earth  model  shown  in  Figure  8.  Since  only  SH-waves  are  radiated  in  this  coordinate 
plane  by  the  source  vector,  the  radial  and  vertical  component  synthetic  seismograms 
contain  no  signal.  Main  reflections  are  identified  as  follows  (prime  indicates  free  surface 
multiple):  (A)  S  to  S  reflection  from  interface  1.  (B)  S  to  S  reflection  from  interface  2. 
(C)  S  to  S  reflection  from  interface  3. 
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Figure  10:  Comparison  of  the  SH-wave  transverse  component  background  displacement  field 
and  the  total  displacement  field  for  the  model  containing  an  isotropic  fracture  zone.  This 
figure  contains  the  time  window  indicated  on  Figure  9.  (A)  Background  displacement 
field.  (B)  Total  displacement  field. 
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Figure  II:  Comparison  of  the  SH  wave  transverse  component  total  displacement  fields  for 
three  orientations  of  the  vertical  fractures  with  respect  tcj  the  receiver  array  from  the 
model  in  Figure  8.  All  plots  in  this  figure  arc  at  the  same  scale  as  those  in  Figure  10.  and 
the  time  window  is  indicated  in  Figure  9.  (A)  Fractures  perpendicular  to  the  receiver 
array.  (B)  Fractures  parallel  to  the  receiver  array.  (C)  Fractures  at  45°  to  the  receiver 
array. 
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Figure  12:  The  radial  component  synthetic  total  field  synthetic  seismogram  resulting  when 
the  anisotropic  fractures  are  aligned  at  4o°  to  the  receiver  array. 
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Figure  13:  Vertical  (above)  and  radial  (below)  component  synthetic  seismograms  from  an 
explosion  source  for  the  background  model  in  Figure  8.  The  main  reflections  are  identified 
as  follows  (prime  indicates  a  free  surface  multiple);  (A)  P  to  P  reflection  from  interface 
I.  (B)  P  to  P  reflection  from  interface  2.  (C)  P  tw  P  reflection  from  interface  3.  (I))  P  to 
S  reflection  from  interface  1.  (F)  P  to  S  reflection  from  interface  2.  (F)  P  to  P  reflection 
from  interface  3.  (C!)  P  to  S  reflection  at  free  surface,  followed  by  .S  to  .S  reflection  at 
interface  3. 
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Figure  14:  Comparison  of  the  P-wave  vertical  component  background  displacement  field  and 
total  displacement  field  for  the  model  containing  an  isotropic  fracture  zone.  Tliis  figure 
contains  the  time  window  indicated  on  Figure  13.  (A)  Background  displacement  field. 
(B)  Total  displacement  field. 
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Figure  16:  Map  showing  the  location  of  the  Larderello  geothermal  field  in  Italy. 
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Figure  17:  East- west  cross-section  of  the  background  earth  model  used  to  compute  the 
synthetic  seismograms  for  the  Badia  1 A  VSP  data.  Sedimentary  units  are  present  in  both 
layers  over  the  L  horizon,  while  the  rock  below  is  metamorphic.  The  seismic  properties 
of  the  metamorphic  zone  arc  relatively  homogeneous  so  that  it  may  be  represented  by  a 
single  velocity.  Compressional  wave  velocities  and  flensities  are  indicated  in  the  figure, 
and  shear  wave  velocities  were  chosen  so  that  tlx*  V^jV,  ratio  was  \/3. 
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Figure  19A:  Positions  of  geophones  and  the  deviated  Badia  lA  well.  The  geophones  are 
|)lotted  in  all  three  coordinate  planes.  A)  Geophone  positions  for  the  zero  offset  source. 
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Figure  19B:  B)  Geophone  positions  for  the  A  offset  source. 
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Figure  22:  Ray  paths  for  a  fan  of  rays  initially  in  the  north-south  vertical  plane.  .Ml  three 
coordinate  planes  are  plotted,  and  the  LI,  L  and  H  markers  are  shown  in  the  appro¬ 
priate  north-south  cross  section.  When  the  rays  encounter  the  LI  and  L  horizon,  the 
irregular  interfaces  bend  the  ray  paths  so  that  they  have  a  strongly  three-dimensional 
propagation  showing  the  complexity  of  wave  propagation  in  the  Larderello  region.  In  a 
one-dimensional  earth  model,  these  rays  would  have  no  east-west  components.  To  com¬ 
pletely  illuminate  the  H  marker,  similar  fans  of  rays  were  traced  over  all  azimuths  from 
l)oth  source  positions.  A)  Ray  fan  from  the  zero  offset  source.  B)  Ray  fan  frc  m  the  .A 
olfset  .source. 


68 


EAST-WEST  (KM) 

0.8  1.0  1.2  1.4  1.6 


EAST-WEST  (KM) 

OJ  1.0  12  1.4  1.6 


EAST-WEST  (KM) 

0.8  1.0  12  1.4  1.6 


EAST-WEST  (KM) 

0.8  1.0  12  1.4  1.6 


Figuri  23:  Contour  plots  of  total  travel  time  from  the  zero  offset  and  A  offset  sources  to 
the  H2  marker  and  then  to  receivers  at  comparable  depths.  Each  contour  line  bounds 
a  Fresnel  zone  as  described  in  the  text.  (A)  Total  travel  time  for  the  zero  offset  source, 
receiver  30.  (B)  Total  travel  time  for  the  zero  offset  source,  receiver  70.  (C)  Total  travel 
lime  for  the  A  offset  source,  receiver  22.  This  receiver  is  at  the  same  depth  as  receiver 
30  for  the  zero  offset  source.  (D)  Total  travel  time  for  the  .A  offset  source  and  receiver 
62.  This  receiver  is  at  the  same  depth  as  receiver  70  for  the  zero  offset  source. 
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